close all
clear all
addpath(genpath(pwd));

%%% Section 1
%%%   Specify important parameters in this section

%% Section 1a
% Dictate folder prefix to save results and figures (will be appended with spreadsheet name):
ImportPlotMasterFolder = 'OUTPUT RevPet OUTPUTS';
SubfolderPrefix = 'RevPet';

% Dictate string to ID model run
SavedName = '';


%% Section 1b INPUT REVPET PARAMETERS
% These parameters are often changed, there are more below in Section 2a
% (search for "Parameters you may want to change (Section 2a)")
% Range of primary melt Mg#s:
MinMgNumPrimary = 0.69;
MaxMgNumPrimary = 0.76;

AcceptableTotalRange = [98.5 101.5];
AcceptableTotalRange = [96 110];

%% Section 1c 
% YOU HAVE TWO CHOICES:
% New Import (A),
% Load .mat file (B),


%if A) IMPORT, then dictate which data should be loaded:

    %if A) IMPORT, then dictate which data should be loaded:
    
    
    %%% Code used for Ontong Java (uncomment to use)
    
%     SpreadsheetName= 'SSS_OntongJava';% Name of spreadsheet
%     TabName = {'PRIMELT3' 'Fitton2004'}; %Spreadsheet Tab(s) with data to plot. DATA IS READ AND PLOTTED IN THIS ORDER
%     DesiredDataName = {'plotOn'}; % keyword(s) that must appear on each row with desired data to plot. If the DesiredDataName is the same for all tabs, you can just specify the one string: e.g.,

    
    %%% Code used to aggregate the 13,539 MORB glasses (uncomment to use)
 
%     SpreadsheetName= 'SSS_13539MORB'; % Name of spreadsheet
%     TabName = {'Gale2013Glasses' 'Marion' 'SGalapagos' 'SGala2' 'S_EPR_PAR' 'MiddleMAR' 'MiddleMAR2' 'MAR3' 'SMAR' 'SEIR' 'CAR' 'AAR'}; %Spreadsheet Tab(s) with data to plot. DATA IS READ AND PLOTTED IN THIS ORDER
%     DesiredDataName = {'REALGALEGAP'}; % keyword(s) that must appear on each row with desired data to plot. If the DesiredDataName is the same for all tabs, you can just specify the one string: e.g.,

    
    %%% Code used to for RevPet Benchmark Transform (uncomment to use)
    
%     SpreadsheetName= 'SSS_Benchmark'; % Name of spreadsheet
%     TabName = {'KnownLherz'}; %Spreadsheet Tab(s) with data to plot. DATA IS READ AND PLOTTED IN THIS ORDER
%     DesiredDataName = {'KnownLherz'}; % keyword(s) that must appear on each row with desired data to plot. If the DesiredDataName is the same for all tabs, you can just specify the one string: e.g.,


    %%% Code used to for Siqueros Transform (comment to use other data)
    
    SpreadsheetName= 'SSS_Siqueros'; % Name of spreadsheet
    TabName = {'Perfit1996' 'Gale2014' 'Gale2013Glasses' 'Gale2013TraceIsotope' 'Niu2008' 'PRIMELT3'}; %Spreadsheet Tab(s) with data to plot:DATA IS READ AND PLOTTED IN THIS ORDER
    DesiredDataName = {'EPR_Compare'}; % keyword(s) that must appear on each row with desired data to plot. If the DesiredDataName is the same for all tabs, you can just specify the one string: e.g.,

    
%% Section 1d
    %if B) LOAD, then dictate the name of the .mat file with data
    Saved_worksheetName2Save = 'SSS_Siqueros.mat';

    % Specify details of new legend, if you are using (can change later in code):
    NewLegendTab = 'Legend_New';
    NewLegendTab_Dataname = 'ReimportLegend'; 
       
       
%%
%%% PROMPTED DECISIONS
%%%   Don't need to change anything here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
addpath(genpath(pwd));



promptreponse123 = input('Do you want to:\n[1] Import from EXCEL \n[2]  Load data from .mat file? \n 1 or 2 [1]:');

if promptreponse123==1
    
    CHOICE = 'A';
    disp('Okay, will import from EXCEL')
    
    clearvars -except SavedName CHOICE subfolder SpreadsheetName TabName DesiredDataName Audio promptreponse123 NewLegendTab NewLegendTab_Dataname  SubfolderPrefix ImportPlotMasterFolder MinMgNumPrimary MaxMgNumPrimary AcceptableTotalRange
    
    
    
    promptreponse = input(sprintf('Is ''%s'' the right spreadsheet to read? \n Press any key if yes, enter ''N'' to change:',SpreadsheetName),'s');
    if strcmp(promptreponse,'N')
        disp('<strong>Ending program. Go manually change variable ''SpreadsheetName'' in MATLAB code</strong>')
        return
    else
        disp('<strong>Great!</strong>')
        disp('     ')
    end
    
    
    
    TabName
    promptreponse = input(sprintf('Are the Excel tab names above correct? \n Press any key if yes, enter ''N'' to change:'),'s');
    if strcmp(promptreponse,'N')
        disp('<strong>Ending program. Go manually change variable ''TabName'' in MATLAB code</strong>')
        return
    else
        disp('<strong>Great!</strong>')
        disp('     ')
    end
    
    
    
    
    if size(DesiredDataName,2)==1
        DesiredDataName_temp=DesiredDataName{1};
        keywordtype = '***All KEYWORDS are the same in each EXCEL tab***';
        DesiredDataName   = TabName;
        for j=1:size(TabName,2)
            DesiredDataName{j} = DesiredDataName_temp;
        end
    else
        keywordtype = '***KEYWORDS are different in each tab***';
    end
    disp('     ')
    disp('     ')
    
    if abs(size(DesiredDataName,2) - size(DesiredDataName,2))>0
        disp(sprintf('<strong>\n Number of EXCEL tabs and datanames do not match. \n Either specify 1 dataname used in all tabs, or specify a unique keyword for each tab. Ending Program</strong>'))
        return
    end
    
    
    
    disp('     ')
    DesiredDataName
    disp(keywordtype)
    disp('     ')
    promptreponse = input(sprintf('Are the above keywords (i.e., datanames) \n that flag the desirable data in each tab correct? \n Press any key if yes, enter ''N'' to change:'),'s');
    if strcmp(promptreponse,'N')
        disp('<strong>Ending program. Go manually change variable ''DesiredDataName'' in MATLAB code</strong>')
        return
    else
        disp('<strong>Great!</strong>')
    end
    
    
    
    
    prompthere = sprintf('\n Will only import major element data with anhydrous total \n between %g%% and %g%% \n KEEP THIS DEFAULT RANGE? \n Press any key if yes, enter ''N'' to change:',...
        AcceptableTotalRange(1),AcceptableTotalRange(2)) ;
        changeRange = input(prompthere,'s');
    
    if strcmp(changeRange,'N')
        AcceptableTotalRange(1) = input('Enter new mininum percentage (e.g., 98.5):');
        AcceptableTotalRange(2) = input('Enter new maxium percentage (e.g., 101.5):');
        fprintf('<strong>Okay, changed to anhydrous totals between %g%% and %g%%</strong>\n',AcceptableTotalRange(1),AcceptableTotalRange(2))
        
    else
        disp('<strong>Okay! Will use default range</strong>')
        
    end
    
    
else if promptreponse123==2
        CHOICE = 'B';
        disp('Okay, will load a previously imported .mat file.')
        
        clearvars -except SavedName CHOICE subfolder SpreadsheetName TabName DesiredDataName Saved_worksheetName2Save Subfolder_previous promptreponse123 ImportPlotMasterFolder SubfolderPrefix MinMgNumPrimary MaxMgNumPrimary
        
        
        
        
        promptreponse = input(sprintf('\n Is ''%s'' the right .mat file to load? \n Press any key if yes, enter ''N'' to change:',Saved_worksheetName2Save),'s');
        if strcmp(promptreponse,'N')
            disp('<strong>Ending program. Go manually change variable ''Saved_worksheetName2Save'' in MATLAB code</strong>')
            return
        else
            disp('<strong>Great! Now importing.</strong>')
            disp('      ')
        end
        
        
    else
        disp('<strong>Ending program. Try again</strong>')
        return
    end
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


if strcmp(CHOICE,'A')==1
    disp('<strong>Now importing...</strong>')
    
    endprogram=0;
    AAA2_Step1_ImportData
    if endprogram ==1
        return
    end
    
    
    
    if size(TabName,2)==1
        worksheetName2Save=sprintf('%s_%s',SpreadsheetName,TabName{:});
    else
        worksheetName2Save=sprintf('%s',SpreadsheetName);
    end
    
    subfolder = sprintf('%s_%s',SubfolderPrefix,worksheetName2Save);
    
    
    promptreponse = input(sprintf('\n Is ''%s'' the right filename to save imported data? \n Press any key if yes, enter ''N'' to change:',worksheetName2Save),'s');
    if strcmp(promptreponse,'N')
        worksheetName2Save = input('Okay, what should it be?','s');
    else
        disp('<strong>Okay! Using that filename</strong>')
    end
    
    

    mkdir([pwd '/' ImportPlotMasterFolder], subfolder);


    
    
    savedetails = input('Do you want to save all (default), or save some (just tabs)?: [all/some]','s');
    if strcmp(savedetails,'some')
        AllVariables = who;
        %find the indicies of a tidbit:
        indicies2keep=[];
        for ek = 1:size(TabName,2)
            IndexC = (strfind(AllVariables,TabName{ek}));
            Index = find((cellfun('isempty', IndexC)==0));
            indicies2keep=[indicies2keep;  Index];
        end
        Variables2save = AllVariables(indicies2keep);
        subfolder='ImportedData';
        eval(['save(''',[pwd '/' ImportPlotMasterFolder '/' subfolder '/' worksheetName2Save '.mat'],''',Variables2save{:})'])
        
    else
        eval(['save(''',[pwd '/' ImportPlotMasterFolder '/' subfolder '/' worksheetName2Save '.mat'],''')'])
    end
    
    
    
    
    
else if strcmp(CHOICE,'B')==1
        disp('<strong>Now loading...</strong>')
        
        SubfolderPrefix_RevPet =SubfolderPrefix;
        ImportPlotMasterFolder_RevPet = ImportPlotMasterFolder;
        
        load(Saved_worksheetName2Save)
        SubfolderPrefix = SubfolderPrefix_RevPet;
        ImportPlotMasterFolder=ImportPlotMasterFolder_RevPet;
        save(Saved_worksheetName2Save,'SubfolderPrefix','ImportPlotMasterFolder','-append')
        
    end
end




if size(TabName,2)==1
    worksheetName2Save=sprintf('%s_%s',SpreadsheetName,TabName{:});
else
    worksheetName2Save=sprintf('%s',SpreadsheetName);
end

if size(DesiredDataName,2)==1
    worksheetName2Save=sprintf('%s_%s',worksheetName2Save,DesiredDataName{1});
end

subfolder = sprintf('%s_%s',SubfolderPrefix,worksheetName2Save);


promptreponse = input(sprintf('\n Is ''%s'' the right filename to save imported data? \n Press any key if yes, enter ''N'' to change:',worksheetName2Save),'s');
if strcmp(promptreponse,'N')
    worksheetName2Save = input('Okay, what should it be?','s');
else
    disp('<strong>Okay! Using that filename</strong>')
end



promptreponse = input(sprintf('\n Is ''%s'' the right subfolder in folder \n ''%s'' to save \n data, reformatted spreadsheet, plots? \n Press any key if yes, enter ''N'' to change:',subfolder,ImportPlotMasterFolder),'s');
if strcmp(promptreponse,'N')
    subfolder = input('Okay, what should it be?','s');
else
    disp('<strong>Okay! Using that folder</strong>')
end

disp('     ')
mkdir([pwd '/' ImportPlotMasterFolder], subfolder);





warning('off','MATLAB:MKDIR:DirectoryExists')


disp('      ')
disp('      ')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Miscll code

%runningColors_FILL = regexprep(runningColors_FILL,'NaN','''k''');
%runningMarkers = regexprep(runningMarkers,'NaN','.');
if promptreponse123==3 || promptreponse123==1
    load('Galeetal2013_SuppTables')
    load('PUM') % Loads trace element source, if you want to change it, change it here!
    ColorsHERE
end


%% Section 2: Set reverse FC path parameters
disp('...')
disp('...')
disp('...')
disp('Now check RevPet settings')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters you *may* want to change (Setion 2a)

MinMgNumOP = 0.56;
MaxMgNumOP = 0.70;

MinMgNumOPA = 0.52;
MaxMgNumOPA = 0.65;

ErOPC= 0.015;
ErOP= 0.015;
PercentAboveCPXsaturation=0.001; %increase to lower threshold for cpx saturation

% ErOPC= 0.02;
% ErOP= 0.02;
ErPrimary = 0.015;

% Dictate name of .mat file with trace element partition coefficients. Can import new ones.
TraceElementPartitionCoeffs = 'Ds_FracCryst45_HiNiD';

% forcing
% MaxMgNumOPA = 0.59;
% MaxMgNumOP = MinMgNumPrimary; %0.70;
% MinMgNumOP = 0.67;
% MinMgNumOPA = 0.57;

prompthere = sprintf('Keep primary melt Mg# range of %g to %g \n Press any key if yes, enter ''N'' to change:',...
    100.*MinMgNumPrimary,100.*MaxMgNumPrimary) ;
changeRange = input(prompthere,'s');
if strcmp(changeRange,'N')
    MinMgNumPrimary = input('Enter new mininum Mg# Primary (e.g., 69):')./100;
    MaxMgNumPrimary = input('Enter new maximum Mg# Primary (e.g., 76):')./100;
    fprintf('<strong>Okay, changed primary melt Mg# range to %g to %g</strong>\n',100.*MinMgNumPrimary,100.*MaxMgNumPrimary)
else
    disp('<strong>Okay! Will use default range</strong>')
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Other parameters you are unlikely to change, but could:

%Specifies which regression to use ("NEW" is Krein et al. 2021,  but there are others including Till et al. 2021, which is "TILL2012"):

whichRegression = 'NEW';
%whichRegression = 'TILL2012';
%whichRegression = 'BBG2020';
%whichRegression = 'Combined';
%whichRegression = 'NEW_noOPX';
%whichRegression = 'Combined_v2';
%whichRegression = 'Combined_Sep';


% Sets how mineral compositions are weighted. Default is even weights
Weights = [1 0 0 0 0];
%Weights = [.4 .15 .15 .2 .15];
% Weights = [0.2 0.2 0.2 0.2 0.2];
% Weights = [1 0 0 0 0];
%     {'Qtz'}    {'Plag'}    {'Oliv'}
%     {'Qtz'}    {'Plag'}    {'Cpx'} dont really use
%     {'Qtz'}    {'Oliv'}    {'Cpx'}
%     {'Plag'}    {'Oliv'}    {'Cpx'}

%vector of pressures of fractional crystallization to test
FracCrystal_P_Try =  [0.001 0.01 0.1:.1:10];

%Don't worry about these at all
distanceValue=2; %this uses OP-OPAM, set to 4 for OP-OPAL   M
percentIncrement = 0.1;
SaturationIncrement = 0.01;


%To force the crystallization of magnetite (or some other mineral, can dictate composition) then use
%this:
MagnetiteInfo = [0.01 .07]; %[MgO content, phase proportion]
Components4OPAM = [1 2 3 4];  %[ 1-Quartz 2-Plag 3-Olivine 4-Cpx]


% For debugging or looking at the optimization of specific points
Figs_OPALM_RMSD = 'false';
SaveFigs_OPALM_RMSD = 'false';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



add2NAME= sprintf('%s_%s_Mg%.3g%.3g_ER%.1g%.1g_%.0f',...
    subfolder,whichRegression, 100*MinMgNumPrimary,100*MaxMgNumPrimary,100*ErOPC,100*ErOP, 100*ErPrimary );
add2NAME= regexprep(add2NAME,'\.','');
worksheetName2Save = [add2NAME];

% Unlikely to change
% initialize and configure name of saved RevPet model

MASTER_REGRESSIONS
clear correctedLiquids FCpath
clear BestFit_FCPath_n test
load(TraceElementPartitionCoeffs);


promptreponse = input(sprintf('Is ''%s'' the right .mat file with trace element partition coeffs? \n Press any key if yes, enter ''N'' to change:',TraceElementPartitionCoeffs),'s');
if strcmp(promptreponse,'N')
    disp('<strong>Ending program. Go manually change variable ''TraceElementPartitionCoeffs'' in MATLAB code</strong>')
    return
    
else
    disp('<strong>Great!</strong>')
end




promptreponse = input(sprintf('\n Do you wish to make an Mg# v. OPLAM diagram for each input basalt? \n Press any key if no, enter ''Y'' to change:'),'s');
if strcmp(promptreponse,'Y')
    Figs_OPALM_RMSD = 'yes';
else
    disp('<strong>Okay! Using that folder</strong>')
end


disp('Everything else is probably okay, but feel free to change other parameters in the code by going to Section 2')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%END OF PROMTED DECISIONS and RevPet Input Settings

%%%%%%%
%% Find optimal FC path
%profile on
close all

disp('Okay, RevPet is calculating...')
data = runningMajors;
DataTrace = runningTrace;
subdata = [1:size(data,1)];
'# of samples'
size(data,1)


%     '# of samples with bad totals'
%     size(find(isnan(subdata(:,3))),1)


% Data imported has columns:
%1-Temp 2-Pressure
% 3-SiO2	4 -TiO2	  5- Al2O3  6-Cr2O3	  7-FeO	  8-MnO	  9-MgO	  10-CaO	11-Na2O   12-K2O
% 13-P2O5 14-NiO 15-H2O 16-total
% 17-Mg# 18-NaK 19 - CaO/Al2O3 20 -
% CaO/Al2O3 moles 21: Na2O/FeO 22: K2O/TiO2 23:1-Mg#
% 24-Quartz 25-Plag 26-Olivine 27-Cpx 28-Oxides 29-Or 30-Ap
Elements = {'T(C)' 'P(kbars)' ...
    'SiO2'	'TiO2'	  'Al2O3' 'Cr2O3'	  'FeO'	  'MnO'	  'MgO'	  'CaO' 'Na2O'   'K2O' 'P2O5' 'NiO' 'H2O' 'total'...
    'Mg#' 'NaK#' 'CaO/Al2O3 wt%' 'CaO/Al2O3 moles' 'Na2O/FeO' 'K2O/TiO2' '1-Mg#' 'Qtz' 'Plag' 'Oliv' 'Cpx' 'Ox' 'Or' 'Ap'};


[A,MgNumCol] = ismember('Mg#', Elements);
[A,NaKNumCol] = ismember('NaK#', Elements);


VariablesT = eval(sprintf('VariablesT_%s',whichRegression));
VariablesX = eval(sprintf('VariablesX_%s',whichRegression));
VariablesP = eval(sprintf('VariablesP_%s',whichRegression));


[A,TemperatureIndicies] = ismember(VariablesT, Elements);
[A,CompIndicies] = ismember(VariablesX, Elements);
[A,PressureIndicies] = ismember(VariablesP, Elements);


FittedTargetVariables = {'Qtz' 'Plag' 'Oliv' 'Cpx'};
[A,FittedVariablesIndicies] = ismember(FittedTargetVariables, Elements);


ModelVariables = {'T' 'Oliv' 'Cpx' 'Plag' 'Qtz'};
[A,ModelVariablesIndicies] = ismember(FittedTargetVariables, ModelVariables);
%TemperatureIndicies, CompIndicies, PressureIndicies, FittedVariablesIndicies, ModelVariablesIndicies

% For Faster Tormey
targetElements ={ 'SiO2'	'TiO2' 'Al2O3'	'Cr2O3'  'FeO'	'MnO'	'MgO'	'CaO'    'Na2O'	'K2O' 'P2O5' 'Fe2O3'};
[A,ElementIndicies4Target_Tormey] = ismember(targetElements, Elements);


FittedTargetVariables_OPALM = {'P' 'T' 'Qtz' 'Plag' 'Oliv' 'Cpx'};
ModelVariables_OPALM = {'P' 'T' 'Oliv' 'Cpx' 'Plag' 'Qtz' 'P'};
[A,ModelVariablesIndiciesOPLAM] = ismember(FittedTargetVariables_OPALM, ModelVariables_OPALM);


targetElements_Yang = {'SiO2' 'TiO2' 'Al2O3' 'Fe2O3' 'FeO' 'MgO' 'CaO' 'Na2O' 'K2O' 'P2O5' 'Cr2O3' 'MnO'};
[A,ElementIndicies4Target_Yang] = ismember(targetElements_Yang, Elements);

ElementsYANG = {'SiO2' 'TiO2' 'Al2O3' 'Fe2O3' 'FeO' 'MgO' 'CaO' 'Na2O' 'K2O' 'P2O5' 'Cr2O3'};
[A,ElementIndicies4Target_OUT] = ismember(Elements, ElementsYANG);


%%

% if interrupted
%  load('Gale_Glass_Trace_24July2019_NMARGlass_Trace_originalKD31')
%  subdata=[401:420];


%% The dual consensus / maximum likelihood algorythm starts here
tic
i=1;
numdata = 0;
%subdata=[9]; %XX
FailMatrix={};
for n=subdata
    
    
    n
    
    
    numdata=numdata+1;
    FailMatrix{numdata}='';
    
    if isnan(data(n,3))==1
        
        % code enters this loop if the major element total is out of bounds
        FailMatrix{numdata}='Bad Total';
        
        GoodEggs_All{numdata,:,:}= NaN.*ones(1,11);
        BagEggs_All{numdata,:,:}=NaN.*ones(1,3);
        
        GoodEggs_ReallyGood_All_Average(numdata,:)=NaN.*zeros(1,12);
        
        SavedResults_v2(numdata,:)=NaN.*zeros(1,10);
        SavedResults_v2_NaN(numdata,:)=NaN.*zeros(1,10);
        fits_optimal(numdata,:)=NaN.*zeros(1,5);
        AllFieldBestFit(numdata,:)=NaN.*zeros(1,3);
        Transitions_v2(numdata,:)=NaN.*zeros(1,2);
        withinPhaseT(numdata,:)=NaN.*zeros(1,3);
        
        BestFitPrimary_n(numdata,:)=NaN.*zeros(1,30);
        BestFitPrimaryTormey_n(numdata,:)=NaN.*zeros(1,7);
        
        BestFitPrimary_Trace_n(numdata,:)=NaN.*zeros(1,45);
        BestFitPrimary_TraceRatios_n(numdata,:)=NaN.*zeros(1,22);
        BestFit_FCPath_Tormey_n(numdata,:)={NaN.*zeros(1,7)};
        BestFit_FCPath_n(numdata,:)={NaN.*zeros(1,30)};
        ReverseFCPath_Trace_n(numdata,:)={NaN.*zeros(1,45)};
        
        Garnet2SpinelTransition(numdata,:) = NaN; %(1.2-0.52.*BestFitPrimary_n(numdata, NaKNumCol)+2.11.*BestFitPrimary_n(numdata,MgNumCol))*10;
        Spinel2PlagTransition(numdata,:) = NaN; %(0.823+1.139.*BestFitPrimary_n(numdata, NaKNumCol)+0.165.*BestFitPrimary_n(numdata,MgNumCol))*10;
        
        bestField_numdata_ALL{numdata,:,:}=NaN.*zeros(1,11);
        bestField_numdata_threshold{numdata,:,:}=NaN.*zeros(1,11);
        
        
        BestFitP(numdata)=NaN;
        OPAM_P_fitness_RMSD(numdata)=NaN;
        PO_OPAM_fitness_RMSD(numdata)=NaN;
        
    else
        clear CorrectedLiquids_ff CorrectedLiquids_zz
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %%% Box "b" in Figure 4
        
        %%% Steps here determine initial viable fractional crystallization paths.
        
        %Save the min Mg# to prevent phases from saturating
        MinMgNum_Data_Real = data(n,MgNumCol);
        MinMgNum_Data = data(n,MgNumCol);
        MinMgNum_Data_Real = MinMgNum_Data;
        MinMgNum_Data=floor(MinMgNum_Data*100)/100;
        
        % the lowest possible Mg# for CPX saturation is the larger of either the MinMgNumOPA as initialized, or the lowest Mg# in the data
        minCpx_Use = max([MinMgNum_Data, MinMgNumOPA]);
        
        % Initial guess of Mg#OPA
        X_targetCPXFoContent_all=[minCpx_Use:SaturationIncrement:MaxMgNumOPA];
        if isempty(X_targetCPXFoContent_all)==1
            X_targetCPXFoContent_all = 0; % If none are viable (if minCpx_Use is >=MaxMgNumOPA), set to zero to avoid OPA
        end
        
        
        % Initial guess of Mg#OP
        minPlag_Use = max([MinMgNum_Data, MinMgNumOP]);
        X_targetPlagFoContent_all = [minPlag_Use:SaturationIncrement:MaxMgNumOP];
        if isempty( X_targetPlagFoContent_all)==1
            X_targetPlagFoContent_all = MinMgNum_Data;% If none are viable (if minPlag_Use is >=MaxMgNumOP), set to the rounded down MinMgNum_Data
        end
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Now we need to check if OPA crystallization is likely by checking
        % to see if the melt is near OPA saturation:
        
        %Finds P of CPX crystallization if it were to have occured
        InitialLiquid_Tormey= TormeyProjection(data(n,:),Elements,ElementIndicies4Target_Tormey);
        PlagOliv_InitialLiquid = PlagOlivProps_QPOC(InitialLiquid_Tormey);
        
        clear TormeyOPAM_Ptest TormeyOPAM_Ptest_big TormeyOPAM_Ptest_big_temp TormeyOPAM_Ptest_MEcomp OPAM_initialLiquid_RMSD fitness_RMSD position_RMSD
        clear TormeyOPAM_Ptest_Save OPALM OPALM_here OPALM_SAVE_Ptest PO_OPALM OPAM_initialLiquid_RMSD PO_OPAM_initialLiquid_RMSD OPALM_initialLiquid_RMSD PO_OPALM_initialLiquid_RMSD DistanceMatrix_NEW_ff
        
        % Finds the best fitting pressure assuming the liquid is cpx saturated
        for P_n =1:size(FracCrystal_P_Try,2)
            OPAM_Ptest_Majors= OPAM_Yang1996(data(n,:),FracCrystal_P_Try(P_n),Elements,ElementIndicies4Target_Yang,ElementIndicies4Target_OUT);
            
            
            TormeyOPAM_Ptest=TormeyProjection(OPAM_Ptest_Majors,Elements,ElementIndicies4Target_Tormey);
            PlagOliv_OPAM = PlagOlivProps_QPOC(TormeyOPAM_Ptest);
            
            TormeyOPAM_Ptest_Save(P_n,:)=TormeyOPAM_Ptest;
            
            [OPALM] = calculate_OPALM_MASTER(data(n,:),FracCrystal_P_Try(P_n),TemperatureIndicies,CompIndicies,ModelVariablesIndiciesOPLAM, ...
                eval(sprintf('garnet_coefficients_%s',whichRegression)), eval(sprintf('spinel_coefficients_%s',whichRegression)), eval(sprintf('plagioclase_coefficients_%s',whichRegression)));
            
            OPALM_here = OPALM(3:end);
            OPALM_SAVE_Ptest(P_n,:)=OPALM_here;
            
            PO_OPALM = PlagOlivProps_QPOC(OPALM_here);
            
            OPAM_initialLiquid_RMSD(P_n,:) = RMSD(InitialLiquid_Tormey,TormeyOPAM_Ptest,Components4OPAM,'NormalizeComponents');
            PO_OPAM_initialLiquid_RMSD(P_n,:) = RMSD(PlagOliv_InitialLiquid,PlagOliv_OPAM,[1 2],'NormalizeComponents');
            
            OPALM_initialLiquid_RMSD(P_n,:) = RMSD(InitialLiquid_Tormey,OPALM_here,Components4OPAM,'NormalizeComponents');
            PO_OPALM_initialLiquid_RMSD(P_n,:) = RMSD(PlagOliv_InitialLiquid,PO_OPALM,[1 2],'NormalizeComponents');
        end
        
        [fitness_RMSD, position_RMSD] = min(OPAM_initialLiquid_RMSD,[],1);
        
        
        BestFitP(n) =  FracCrystal_P_Try(position_RMSD)   ;
        OPAM = TormeyOPAM_Ptest_Save(position_RMSD,:);
        OPAM_P = BestFitP(n);
        OPAM_P_fitness_RMSD(n) = fitness_RMSD;
        
        
        [PO_fitness_RMSD, PO_position_RMSD] = min(PO_OPAM_initialLiquid_RMSD,[],1);
        PO_OPAM_BestFitP =  FracCrystal_P_Try(PO_position_RMSD);
        PO_OPAM_P = PO_OPAM_BestFitP;
        PO_OPAM_fitness_RMSD(n) = PO_fitness_RMSD;
        
        
        [OPALM_fitness_RMSD, OPALM_position_RMSD] = min(OPALM_initialLiquid_RMSD,[],1);
        OPALM_BestFitP =  FracCrystal_P_Try(OPALM_position_RMSD);
        OPALM_P = OPALM_BestFitP;
        OPALM_P_fitness_RMSD = OPALM_fitness_RMSD;
        
        [PO_OPALM_fitness_RMSD, PO_OPALM_position_RMSD] = min(PO_OPALM_initialLiquid_RMSD,[],1);
        PO_OPALM_BestFitP =  FracCrystal_P_Try(PO_OPALM_position_RMSD);
        PO_OPALM_P = PO_OPALM_BestFitP;
        PO_OPALM_P_fitness_RMSD = PO_OPALM_fitness_RMSD;
        
        
        OP_OPAM = PlagOlivProps_QPOC(OPAM);
        OP_OPAM_ALL_P =  PlagOlivProps_QPOC(TormeyOPAM_Ptest_Save);
        OP_OPALM_ALL_P =  PlagOlivProps_QPOC(OPALM_SAVE_Ptest);
        
        
        %'Determine if liquid could be saturated in OP'
        % CHECK OPALM point:
        a_PO_OPALM=find(OP_OPALM_ALL_P(:,1)<PlagOliv_InitialLiquid(1) & PO_OPALM_initialLiquid_RMSD < ErOP);
        a_PO_OPAM=find(PO_OPAM_initialLiquid_RMSD < ErOP);
        
        % Determine if liquid could be saturated in OPC
        % Needs to be close to OPAM point AND sufficiently saturated in CPX
        % PercentAboveCPXsaturation: increase PercentAboveCPXsaturation to lower threshold for cpx saturation
        if round(InitialLiquid_Tormey(4),4) >= round(OPAM(4)*(1-PercentAboveCPXsaturation),4)
            a_OPAM = find(OPAM_initialLiquid_RMSD < ErOPC);
            if any(a_OPAM) ==1
                Pressures_OPA =  [min(FracCrystal_P_Try(a_OPAM)) max(FracCrystal_P_Try(a_OPAM))];
                Pressures_OP = Pressures_OPA;
            else
                X_targetCPXFoContent_all = 0;
                Pressures_OP =  [min(FracCrystal_P_Try(a_PO_OPAM)) max(FracCrystal_P_Try(a_PO_OPAM))];
                if any(Pressures_OP)==0
                    Pressures_OP=[0 0];
                end
                
            end
            
        else
            X_targetCPXFoContent_all = 0;
            Pressures_OP =  [min(FracCrystal_P_Try(a_PO_OPAM)) max(FracCrystal_P_Try(a_PO_OPAM))];
            if any(Pressures_OP)==0
                Pressures_OP=[0 0];
            end
        end
        
        
        clear MasterSaturationMgNums
        xxx=1;
        for qw = 1:size(X_targetCPXFoContent_all,2)
            for qw2 = 1:size(X_targetPlagFoContent_all,2)
                MasterSaturationMgNums(xxx,1) = X_targetCPXFoContent_all(qw);
                MasterSaturationMgNums(xxx,2) = X_targetPlagFoContent_all(qw2);
                xxx=xxx+1;
            end
        end
        
        % remove fractionation paths that have cpx coming in before plag.
        % Future work can easily implement this functionality in ! (But please check and calibrate!)
        cpxBeforePlag = find(MasterSaturationMgNums(:,2)<MasterSaturationMgNums(:,1));
        MasterSaturationMgNums(cpxBeforePlag,:)=[];
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %%% Box "c" in Figure 4
        
        clear ReverseFCPath_zz
        clear SpinelFits
        clear PrimitveMelts_fitsPTTmp_zz
        clear FCparameters_zz
        clear PrimitveMeltsBestFit_zz
        clear PrimitveMeltsFits_zz
        clear TormeyPrimary_zz
        clear PrimitveMelts
        clear CorrectedLiquids_zz
        clear BestFit2Optimize
        clear correctedLiquids
        clear Data_all
        clear Data
        clear TormeyData
        clear results
        clear fits
        clear fitsPTTmp_output
        clear GoodEggs BadEggs    ReverseFCPath_zz Fout_save phaseChangeInfo_best
        
        GoodEggs=NaN.*ones(1,11);
        BadEggs = NaN.*ones(1,3);
        
        
        zz=1;
        badCount = 1;
        
        % If you wanted to manually change the tested fractionation paths,
        % redefine it here by uncommenting the example below:
        %  MasterSaturationMgNums=[.65 .67];
        
        for ff = 1:size(MasterSaturationMgNums,1)
            
            DistanceMatrix_NEW_ff=[];
            
            MgNumOPA = MasterSaturationMgNums(ff,1);
            MgNumOP = MasterSaturationMgNums(ff,2);
            
            %"RFC_function" function actually does the work of calculating a
            %single reverse fractional crystallization path :
            
            [CorrectedLiquids_ff, cpx_ff, temperature_ff,phaseProportions_out_ff, phaseChangeInfo_ff,AllCorrectedLiquids_Trace_ff,...
                solidphases_ff, DistanceMatrix_ff,OPAMout_ff,OPALMout_ff,OlivineKDused_ff,liquidsOP_ff,Fout_ff]...
                =  ...
                RFC_function(data(n,:), OPAM_P,Pressures_OP,...
                percentIncrement,MgNumOPA, MgNumOP, MaxMgNumPrimary,ErOPC,ErOP, ...
                MagnetiteInfo,whichRegression,distanceValue,Elements,Weights,Components4OPAM,...
                TemperatureIndicies, CompIndicies, PressureIndicies, FittedVariablesIndicies, ModelVariablesIndicies,ElementIndicies4Target_Tormey,...
                ModelVariablesIndiciesOPLAM,ElementIndicies4Target_Yang,ElementIndicies4Target_OUT,...
                DataTrace(n,:),Dmatrix);
            
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            %%% Box "d" in Figure 4. These steps find the most likely primary
            %%% melt along the given ff path
            
            
            if (isnan(CorrectedLiquids_ff))==1
                % Entering this loop means the tested Mg#OP and Mg#OPA failed to stay close
                % to required phase boundaries
                BadEggs(badCount,:) = [ff MgNumOPA MgNumOP];
                badCount = badCount+1;
                
            else
                % Enter this loop if there are viable RFPs
                
                
                
                % useful figure showing distances to all saturation boundaries as a function of Mg#
                %                                 figure(300+n)
                %                                 close
                %                                 figure(300+n)
                %                                 hold on
                %                                 plot(CorrectedLiquids_ff(:,MgNumCol), DistanceMatrix_ff(:,1), 'o-k','Color',cadmiumgreen);
                %                                 plot(CorrectedLiquids_ff(:,MgNumCol), DistanceMatrix_ff(:,3), 'o-k','Color',coffee);
                %                                 plot(CorrectedLiquids_ff(:,MgNumCol), DistanceMatrix_ff(:,2), '.:b','color',cadmiumgreen,'LineWidth',4);
                %                                 plot(CorrectedLiquids_ff(:,MgNumCol), DistanceMatrix_ff(:,4), '.:b','color',coffee,'LineWidth',4);
                %                                 plot(CorrectedLiquids_ff(:,MgNumCol), DistanceMatrix_ff(:,5), '-k','markerfacecolor',black,'Color',black,'LineWidth',4);
                %                                 plot(CorrectedLiquids_ff(:,MgNumCol), DistanceMatrix_ff(:,6), '--k','markerfacecolor',black,'Color',black,'LineWidth',4);
                %                                 plot(CorrectedLiquids_ff(:,MgNumCol), DistanceMatrix_ff(:,7), ':k','markerfacecolor',black,'Color',black,'LineWidth',4);
                %                                 yline(ErOPC);
                %                                 yline(ErOP,'b');
                %                                 legend('OPAM (:,1)', 'OPALM (:,3)', 'OPM from OPAM (:,2)', 'OPM from OPALM (:,4)', 'closest melting OPALM garnet','closest melting OPALM spinel','closest melting OPALM plag')
                %                                 set(gcf,'position',[560 528 560 420])
                
                
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                
                
                % these steps tag failures that have pressure predictions outside
                % of phase stability (i.e., garnet peridotite at 8 kbars)
                
                %Save melting temps in a separate matrix
                T_temp_DistanceMatrix_ff=DistanceMatrix_ff(:,[14 15 16]);
                
                % Determine phase transitions within error
                % All solutions with melting pressures outside acceptable
                % ranges will be considered failures:
                
                Gar2SpTrans = DistanceMatrix_ff(:,9)-1;%0.81;
                Sp2PlagTrans_Sp = DistanceMatrix_ff(:,10)-2.5;%1
                Sp2PlagTrans_Plag = DistanceMatrix_ff(:,10)+1;
                
                
                % Calculate melting pressure distances from phase transitions
                DistanceMatrix_ff(:,14) = DistanceMatrix_ff(:,11)-Gar2SpTrans;
                
                DistanceMatrix_ff(:,15) = Gar2SpTrans-DistanceMatrix_ff(:,12);
                DistanceMatrix_ff(:,16) = DistanceMatrix_ff(:,12)-Sp2PlagTrans_Sp;
                
                DistanceMatrix_ff(:,17) = Sp2PlagTrans_Plag-DistanceMatrix_ff(:,13);
                DistanceMatrix_ff = [DistanceMatrix_ff T_temp_DistanceMatrix_ff];
                
                F= Fout_ff(:,2);
                Fo= CorrectedLiquids_ff(:,MgNumCol)*100;
                
                
                % Create a new distance matrix that removes any possible solutions
                % that do not fall within the predicted phase transition boundaries
                % (within the error set above)
                passGar = DistanceMatrix_ff(:,14);
                passGar(passGar<0) = -1;
                passGar(passGar>0) = 1;
                
                
                passSpin = DistanceMatrix_ff(:,15).*DistanceMatrix_ff(:,16);
                passSpin(passSpin<0) = -1;
                passSpin(passSpin>0) = 1;
                
                
                
                passPlag = DistanceMatrix_ff(:,17);
                passPlag(passPlag<0) = -1;
                passPlag(passPlag>0) = 1;
                
                
                %extract the distances from the OPALM points for each field
                DistanceMatrix_NEW_ff(:,1)=DistanceMatrix_ff(:,5).*passGar;
                DistanceMatrix_NEW_ff(:,2)=DistanceMatrix_ff(:,6).*passSpin;
                DistanceMatrix_NEW_ff(:,3)=DistanceMatrix_ff(:,7).*passPlag;
                
                
                % Create new distance matrix that takes into account the phase
                % transition requirement
                DistanceMatrix_NEW_ff_2 = DistanceMatrix_NEW_ff;
                DistanceMatrix_NEW_ff_2(DistanceMatrix_NEW_ff_2<=0)=nan;
                DistanceMatrix_NEW_ff_fo_2 = DistanceMatrix_NEW_ff_2;
                DistanceMatrix_NEW_ff_fo_2(Fo<(100*MinMgNumPrimary),:)=nan;
                
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

                % now actually find the most likely primary melt along the
                % path
                
                % Find which melts have high enough Mg#s to be primary
                FinalFoTest=find(CorrectedLiquids_ff(:,MgNumCol)>MinMgNumPrimary);
                
                %Set any solutions with Mg# lower than MinMgNumPrimary to be
                %failures
                DistanceMatrix_NEW_ff_2(1:FinalFoTest(1)-1,:) = nan;
                
                
                %Save the minimum melting OPLAM distance across all fields
                min_OPALM_MELT_ALL = min(DistanceMatrix_ff(:,5:7),[],2); %This finds at all melts along the path
                min_OPALM_MELT = min(DistanceMatrix_NEW_ff_2,[],2); %This finds only those with primary melt mg#s
                
                %Save the minimum melting OPLAM distance in each field
                %assuming:
                [EachFieldMin_OPALM_ALL,EachFieldMin_OPALM_Index_ALL]=min(DistanceMatrix_NEW_ff_2,[],1);
                
                
                % Find which melts have low enough OPALM error and high enough
                % Mg#s:
                [a_O,b]=find(min_OPALM_MELT < ErPrimary & CorrectedLiquids_ff(:,MgNumCol)>MinMgNumPrimary);
                
                
                
                % BestFit2Optimize_lastFo=min_OPALM_MELT';
                if any(a_O) == 1
                    %Enter this loop if there is one or more RFPs that come
                    %withing ErPri of a MSP
                    
                    %This code then assumes the best fit primary melt is the
                    %one in the group with the lowest Mg#
                    
                    
                    FirstFoIndex = a_O(1);
                    FirstFo = CorrectedLiquids_ff(FirstFoIndex,MgNumCol);
                    
                    MinError = min_OPALM_MELT(FirstFoIndex);
                    %                min_OPALM_MELT2 = min(min_OPALM_MELT);
                    %                [b,c]=min(min_OPALM_MELT2);
                    
                else
                    
                    %Enter this loop if no RFP comes within Er of a MSP
                    %In this case assume that the "best fit" primary melt is the
                    %one in the group with the lowest error (i.e. lowest
                    %RMSD_OPALM)
                    
                    min_OPALM_MELT2 = min(DistanceMatrix_NEW_ff_2(FinalFoTest,:),[],2);
                    [MinError,c]=min(min_OPALM_MELT2);
                    
                    
                    FirstFoIndex = FinalFoTest(c);
                    FirstFo = CorrectedLiquids_ff(FirstFoIndex,MgNumCol);
                end
                
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                
                
                % Work in progress to save the best fitting primary melt
                % for each stablity field for statistics. I moved this
                % around a bit.
                
                
                ErOkay = 0.02;
                DistanceMatrix_NEW_ff_fo_2_Okay= DistanceMatrix_NEW_ff_fo_2;
                DistanceMatrix_NEW_ff_fo_2_Okay(DistanceMatrix_NEW_ff_fo_2_Okay>ErOkay) = NaN;
                %Find the minimum Er in each field
                % https://www.mathworks.com/matlabcentral/answers/69517-find-first-value-below-a-minimum-in-a-vectorized-way
                A = DistanceMatrix_NEW_ff_fo_2_Okay';
                thresh = ErOkay;
                temp_mat = repmat(1:size(A,2),size(A,1),1);
                temp_mat = temp_mat .* (A<=thresh);
                temp_mat(temp_mat == 0) = NaN;
                [your_idx,your_idx2] = min(temp_mat,[],2);
                EachFieldMin_OPALM = your_idx';
                EachFieldMin_OPALM_Index = your_idx2';
                ErField_okay_ff = diag(DistanceMatrix_NEW_ff_fo_2_Okay(EachFieldMin_OPALM_Index,:))';
                ErField_okay_ff(EachFieldMin_OPALM==NaN) = NaN;
                ErField_okay_ff_index = EachFieldMin_OPALM;
                
                %Find the minimum Er at each Fo
                [s,k]=min(DistanceMatrix_NEW_ff_fo_2_Okay,[],2);
                DistanceMatrix_NEW_ff_fo_2_Okay = [DistanceMatrix_NEW_ff_fo_2_Okay k s];
                
                
                DistanceMatrix_NEW_ff_fo_2_ErPrimary= DistanceMatrix_NEW_ff_fo_2;
                DistanceMatrix_NEW_ff_fo_2_ErPrimary(DistanceMatrix_NEW_ff_fo_2_ErPrimary>ErPrimary) = NaN;
                %Find the minimum Er in each field
                % https://www.mathworks.com/matlabcentral/answers/69517-find-first-value-below-a-minimum-in-a-vectorized-way
                A = DistanceMatrix_NEW_ff_fo_2_ErPrimary';
                thresh = ErOkay;
                temp_mat = repmat(1:size(A,2),size(A,1),1);
                temp_mat = temp_mat .* (A<=thresh);
                temp_mat(temp_mat == 0) = NaN;
                [your_idx,your_idx2] = min(temp_mat,[],2);
                EachFieldMin_OPALM = your_idx';
                EachFieldMin_OPALM_Index = your_idx2';
                ErField_best_ff = diag(DistanceMatrix_NEW_ff_fo_2_ErPrimary(EachFieldMin_OPALM_Index,:))';
                ErField_best_ff(EachFieldMin_OPALM==NaN) = NaN;
                ErField_best_ff_index = EachFieldMin_OPALM;
                
                
                
                %Find the minimum Er at each Fo
                [s,k]=min(DistanceMatrix_NEW_ff_fo_2_ErPrimary,[],2);
                DistanceMatrix_NEW_ff_fo_2_ErPrimary = [DistanceMatrix_NEW_ff_fo_2_ErPrimary k s];
                
                
                
                Table_DistanceMatrix_ff = table(Fo, F, DistanceMatrix_ff, passGar,passSpin,passPlag,DistanceMatrix_NEW_ff_2,DistanceMatrix_NEW_ff_fo_2);
                varNames1 = {'ER_OPAM' 'ER_OP_OPAM' 'ER_OPALM' 'ER_OP_OPALM' 'ER_G' 'ER_S' 'ER_P'...
                    'AboveCPX' 'G2S' 'S2P' 'PG' 'PS' 'PP' 'GG2S' 'G2SS' 'SS2P' 'S2PP' 'TG' 'TS' 'TP'};
                varNames2 = {'ER_G_passP' 'ER_S_passP' 'ER_P_passP'};
                varNames3 = {'ER_G_passP_Fo' 'ER_S_passP_Fo' 'ER_P_passP_Fo'};
                Table_DistanceMatrix_ff = splitvars(Table_DistanceMatrix_ff,'DistanceMatrix_ff','NewVariableNames',varNames1);
                Table_DistanceMatrix_ff = splitvars(Table_DistanceMatrix_ff,'DistanceMatrix_NEW_ff_2','NewVariableNames',varNames2);
                Table_DistanceMatrix_ff = splitvars(Table_DistanceMatrix_ff,'DistanceMatrix_NEW_ff_fo_2','NewVariableNames',varNames3);
                
                
                % https://www.mathworks.com/matlabcentral/answers/69517-find-first-value-below-a-minimum-in-a-vectorized-way
                A = DistanceMatrix_NEW_ff_2';
                thresh = ErPrimary;
                temp_mat = repmat(1:size(A,2),size(A,1),1);
                temp_mat = temp_mat .* (A<=thresh);
                temp_mat(temp_mat == 0) = NaN;
                [your_idx,your_idx2] = min(temp_mat,[],2);
                EachFieldMin_OPALM = your_idx';
                EachFieldMin_OPALM_Index = your_idx2';
                
                
                
                %             EachFieldMin_OPALM_Index
                %
                %             fixHighEr = find(EachFieldMin_OPALM_Index==1);
                %             for ssh = fixHighEr
                %                 [EachFieldMin_OPALMaa,EachFieldMin_OPALM_Indexaa]=min(abs(DistanceMatrix_NEW_ff(:,ssh)));
                %                 EachFieldMin_OPALM_Index(ssh)=EachFieldMin_OPALM_Indexaa;
                %             end
                %             EachFieldMin_OPALM_Index
                
                
                
                Er_ThreshGSP =  [DistanceMatrix_NEW_ff_2(EachFieldMin_OPALM_Index(1),1) DistanceMatrix_NEW_ff_2(EachFieldMin_OPALM_Index(2),2) DistanceMatrix_NEW_ff_2(EachFieldMin_OPALM_Index(3),3)];
                Fo_ThreshGSP = CorrectedLiquids_ff(EachFieldMin_OPALM_Index,MgNumCol)';
                P_ThreshGSP = [DistanceMatrix_ff(EachFieldMin_OPALM_Index(1),11), DistanceMatrix_ff(EachFieldMin_OPALM_Index(2),12), DistanceMatrix_ff(EachFieldMin_OPALM_Index(3),13)];
                F_ThresGSP = Fout_ff(EachFieldMin_OPALM_Index,2)';
                
                Table_EachFieldBestFit_zz_thresh(zz,:) = table(ff, MgNumOPA, MgNumOP,Er_ThreshGSP,Fo_ThreshGSP,P_ThreshGSP,F_ThresGSP);
                
                
                EachFieldBestFit_zz_threshold(zz,:) = table2array(Table_EachFieldBestFit_zz_thresh(zz,:) );
                %                 EachFieldBestFit_zz_threshold(zz,:) = [ff MgNumOPA MgNumOP...
                %                   Er_ThreshGSP...
                %                     CorrectedLiquids_ff(EachFieldMin_OPALM_Index,MgNumCol)' ...
                %                     DistanceMatrix_ff(EachFieldMin_OPALM_Index(1),11) DistanceMatrix_ff(EachFieldMin_OPALM_Index(2),12) DistanceMatrix_ff(EachFieldMin_OPALM_Index(3),13)...
                %                     Fout_ff(EachFieldMin_OPALM_Index,2)'];
                
                
                EachFieldBestFit_zz_All(zz,:) = [ff MgNumOPA MgNumOP...
                    DistanceMatrix_NEW_ff(EachFieldMin_OPALM_Index_ALL(1),1) DistanceMatrix_NEW_ff(EachFieldMin_OPALM_Index_ALL(2),2) DistanceMatrix_NEW_ff(EachFieldMin_OPALM_Index_ALL(3),3)...
                    CorrectedLiquids_ff(EachFieldMin_OPALM_Index_ALL,MgNumCol)' ....
                    DistanceMatrix_ff(EachFieldMin_OPALM_Index_ALL(1),11) DistanceMatrix_ff(EachFieldMin_OPALM_Index_ALL(2),12) DistanceMatrix_ff(EachFieldMin_OPALM_Index_ALL(3),13)...
                    Fout_ff(EachFieldMin_OPALM_Index_ALL,2)'];
                
                
                % end work in progress
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                
               
                
                
                
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                
                %%% "Box e in Figure 4"  Save results for this particular RFP!
                
                
                CorrectedLiquids_zz(zz,:)=CorrectedLiquids_ff(FirstFoIndex,:);
                TormeyPrimary_zz(zz,:)= TormeyProjection(CorrectedLiquids_ff(FirstFoIndex,:),Elements,ElementIndicies4Target_Tormey);
                CorrectedLiquids_Trace_zz(zz,:)=AllCorrectedLiquids_Trace_ff(FirstFoIndex,:);
                
                
                
                ReverseFCPath_zz{zz,:,:}=CorrectedLiquids_ff(1:FirstFoIndex,:);
                Fout_zz{zz,:,:}=Fout_ff;
                ReverseFCPath_Trace_zz{zz,:,:}=AllCorrectedLiquids_Trace_ff(1:FirstFoIndex,:);
                
                
                % Rerun XPTmeter_MASTER (the function that calculates
                % compositional distance to the best fitting primary melt
                % to save more results)
                
                %XPTmeter_MASTER, which is also called in the function
                %RFC_function, takes in the melt comp data, and calculates the
                %distance to the nearest OPLAM points using expressions as
                %in Table 1.  The function however can use other
                %regressions, such as those from Grove et al. (2013) and
                %Till et al. (2012).  The regression choice is changed by changing the variable "whichRegression"
                
                results = XPTmeter_MASTER(CorrectedLiquids_zz(zz,:),TemperatureIndicies, CompIndicies, PressureIndicies, FittedVariablesIndicies, ModelVariablesIndicies,...
                    eval(sprintf('garnet_coefficients_%s',whichRegression)), eval(sprintf('spinel_coefficients_%s',whichRegression)), eval(sprintf('plagioclase_coefficients_%s',whichRegression)));
                
                
                if abs(size(DistanceMatrix_NEW_ff_2,1)- size(CorrectedLiquids_ff,1)) >0
                    'YIKES FIX'
                    return
                end
                
                
                fits = DistanceMatrix_NEW_ff_2(FirstFoIndex,:);
                fitsPTTmp_output = results.PTTmp;
                
                %find the best fit composition from any of the fields:
                % [NRMSD_BEST BestFitFIELD]
                
                [fits(:,4),fits(:,5)]=min(fits,[],2);
                BestFit2Optimize = fits(:,4);
                fitField = fits(:,5);
                
                
                PrimitveMeltsFits_zz(zz,:) = fits;
                PrimitveMeltsBestFit_zz(zz) = BestFit2Optimize;
                FCparameters_zz(zz,:)=[MgNumOPA,MgNumOP];
                PrimitveMelts_fitsPTTmp_zz(zz,:) = fitsPTTmp_output;
                
                phaseChangeInfo_zz(zz,:,:) = phaseChangeInfo_ff;
                
                %find the best fit composition from only the spinel field:
                SpinelFits(zz) =  fits(:,2);
                
                
                if fitField == 1
                    BestFitPT =fitsPTTmp_output(1:3);
                else if fitField == 2
                        BestFitPT= fitsPTTmp_output(4:6);
                    else if fitField == 3
                            BestFitPT= fitsPTTmp_output(7:9);
                        end
                    end
                end
                
                % This variable saves the details of any RFP that "works"
                GoodEggs(zz,:) = [ff MgNumOPA MgNumOP FirstFo MinError fitField BestFitPT Fout_ff(FirstFoIndex,2) FirstFoIndex];
                
                
                
                %%%
                %%%
                %%%
                %%% useful plots, examples are in Figure 4.  Each basalt has its own plot.
                if strcmp(Figs_OPALM_RMSD,'yes')==1
                    figure(numdata)
                    symbols4plot = {'kd','kv','ks'};
                    hold on
                    axis square
                    %axis([100*MinMgNumPrimary 100*MaxMgNumPrimary+1 0 8])
                    %axis([100*MinMgNum_Data 100*MaxMgNumPrimary+1 0 8])
                    ytickformat(gca, 'percentage');
                    
                    if zz==1
                        xline(MinMgNum_Data_Real*100,'HandleVisibility','off');
                        yline(ErPrimary*100,'k-','Er_{Primary}','LabelVerticalAlignment','middle','FontName','Times New Roman','FontSize',12);
                    end
                    
                    %plot(100.*CorrectedLiquids_ff(:,MgNumCol),  100.*min_OPALM_MELT,'.-','Color',red)
                    %plot(100.*CorrectedLiquids_ff(:,MgNumCol),  100.*min_OPALM_MELT_ALL,'.-','Color',grey9)
                    
                    plot(100.*CorrectedLiquids_ff(:,MgNumCol),  100.*min_OPALM_MELT_ALL,'.-','Color','r')
                    
                    plot(100.*CorrectedLiquids_ff(1:min(FinalFoTest)-1,MgNumCol),  100.*min_OPALM_MELT_ALL(1:min(FinalFoTest)-1),'.-','Color',grey9)
                    
                    
                    
                    MinMgNum_Data
                    MgNumOPA
                    MgNumOP
                    
                    if MgNumOPA>=MinMgNum_Data_Real
                        plot(100.*CorrectedLiquids_ff(FinalFoTest,MgNumCol),  100.*min_OPALM_MELT(FinalFoTest),'k.-','Color',cadmiumgreen,'linewidth',2,'markersize',15)
                        plot(100.*FirstFo, 100.*MinError,symbols4plot{fits(5)},'MarkerFaceColor',cadmiumgreen,'MarkerSize',13)
                    else if MgNumOP>=MinMgNum_Data_Real
                            plot(100.*CorrectedLiquids_ff(FinalFoTest,MgNumCol),  100.*min_OPALM_MELT(FinalFoTest),'k.-','Color',airsuperiorityblue,'linewidth',2,'markersize',15)
                            plot(100.*FirstFo, 100.*MinError,symbols4plot{fits(5)},'MarkerFaceColor',airsuperiorityblue,'MarkerSize',13)
                        else
                            plot(100.*CorrectedLiquids_ff(FinalFoTest,MgNumCol),  100.*min_OPALM_MELT(FinalFoTest),'k.-','linewidth',2,'markersize',15)
                            plot(100.*FirstFo, 100.*MinError,symbols4plot{fits(5)},'MarkerFaceColor',grey5,'MarkerSize',13)
                        end
                    end
                    
                    
                    
                    %             if FirstFo > MaxMgNumPrimary
                    %                 axis([100*MinMgNumPrimary 100*FirstFo+1 0 8])
                    %             end
                    set(gcf,'position',[177 223 743 669])
                    set(gca,'fontsize', 20,'LineWidth',.7,'FontName','Times New Roman')
                    box on
                    ax=gca;
                    set(gca,'XColor', 'k')
                    set(gca,'YColor', 'k')
                    set(gca, 'Xtick',0:1:100)
                    ax.XAxis.MinorTick = 'on';
                    ax.XAxis.MinorTickValues = ax.XLim(1):0.5:ax.XLim(2);
                    ax.YAxis.MinorTick = 'on';
                    ax.YAxis.MinorTickValues = ax.YLim(1):.5:ax.YLim(2);
                    set(gca,'ticklength',1.5*[0.0200    0.0500])
                    xlabel('Mg# Melt')
                    ylabel('RMSD to nearest OPALM')
                    titleHERE = sprintf('P%s-%s-%s',num2str(n),CombinedLabel{n},whichRegression);
                    set(gcf,'name',titleHERE)
                    titleleg = 'blk=O,blue=OP,grn=OPA;sq=plg,tri=spn,diam=gar';
                    title(titleleg);
                end
                %%%
                %%%
                %%%
                %%%
                %%%
                
                zz=zz+1;
                
            end
            
            
            if mod(ff,50)==0
                ff
            end
            
        end
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        %%% "Box f and g in Figure 4"
        % All possible RFP have now been tested, and we found the most
        % likely pimary melt for each path.
        
        % Now we find which path of all the RFPs is best
        
        % Other ways of saying the same thing:
        % FINISHED SEARCHING ff paths.  Now find of the successful zz paths, which one best fits a multiple saturation point
        % Finds the best fit FC path for each basalt (n)
        
        
        if any(GoodEggs)==0
            % Enter this loop if all RFPs failed.
            FailMatrix{numdata}='No FC Paths';
            
            GoodEggs_All{numdata,:,:}= NaN.*ones(1,11);
            BagEggs_All{numdata,:,:}=NaN.*ones(1,3);
            
            GoodEggs_ReallyGood_All_Average(numdata,:)=NaN.*zeros(1,12);
            
            SavedResults_v2(numdata,:)=NaN.*zeros(1,10);
            SavedResults_v2_NaN(numdata,:)=NaN.*zeros(1,10);
            fits_optimal(numdata,:)=NaN.*zeros(1,5);
            AllFieldBestFit(numdata,:)=NaN.*zeros(1,3);
            Transitions_v2(numdata,:)=NaN.*zeros(1,2);
            withinPhaseT(numdata,:)=NaN.*zeros(1,3);
            
            BestFitPrimary_n(numdata,:)=NaN.*zeros(1,30);
            BestFitPrimaryTormey_n(numdata,:)=NaN.*zeros(1,7);
            
            BestFitPrimary_Trace_n(numdata,:)=NaN.*zeros(1,45);
            BestFitPrimary_TraceRatios_n(numdata,:)=NaN.*zeros(1,22);
            BestFit_FCPath_Tormey_n(numdata,:)={NaN.*zeros(1,7)};
            BestFit_FCPath_n(numdata,:)={NaN.*zeros(1,30)};
            ReverseFCPath_Trace_n(numdata,:)={NaN.*zeros(1,45)};
            
            Garnet2SpinelTransition(numdata,:) = NaN; %(1.2-0.52.*BestFitPrimary_n(numdata, NaKNumCol)+2.11.*BestFitPrimary_n(numdata,MgNumCol))*10;
            Spinel2PlagTransition(numdata,:) = NaN; %(0.823+1.139.*BestFitPrimary_n(numdata, NaKNumCol)+0.165.*BestFitPrimary_n(numdata,MgNumCol))*10;
            
            bestField_numdata_ALL{numdata,:,:}=NaN.*zeros(1,11);
            bestField_numdata_threshold{numdata,:,:}=NaN.*zeros(1,11);
            
            BestFitP(numdata)=NaN;
            OPAM_P_fitness_RMSD(numdata)=NaN;
            PO_OPAM_fitness_RMSD(numdata)=NaN;
            
        else
            
            
            
            [optimalbestfit,optimalbestfit_index] = min(PrimitveMeltsBestFit_zz);
            
            % If errors are really low, pick the first Mg# that comes within error
            if optimalbestfit < ErPrimary
                optimalbestfit_index=find(PrimitveMeltsBestFit_zz<ErPrimary,1,'first');
                optimalbestfit=    PrimitveMeltsBestFit_zz(optimalbestfit_index);
                
            end
            
            [optimalbestfit_SPINEL,optimalbestfit_index_SPINEL]= min(SpinelFits);
            
            phaseChangeInfo_best = squeeze(phaseChangeInfo_zz(optimalbestfit_index,:,:));
            Fout_best = squeeze(Fout_zz{optimalbestfit_index,:,:});
            CorrectedLiquids_ff = squeeze(ReverseFCPath_zz{optimalbestfit_index,:,:});
            
            
            phaseChangeInfo_best(3) = Fout_best(GoodEggs(optimalbestfit_index,11),2) - phaseChangeInfo_best(3);
            phaseChangeInfo_best(4) = GoodEggs(optimalbestfit_index,10);
            
            phaseChangeInfo_best(phaseChangeInfo_best<0)=0; % XX SHOULD FIX!!!
            
            phaseChangeInfo_n(numdata,:)  = phaseChangeInfo_best;
            
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            %Saves the best FC path
            BestFit_FCPath_n{numdata,:,:}=ReverseFCPath_zz{optimalbestfit_index,:,:};
            ReverseFCPath_Trace_n{numdata,:,:}=ReverseFCPath_Trace_zz{optimalbestfit_index,:,:};
            
            OptimalFCPath_SPINEL{numdata,:,:}=ReverseFCPath_zz{optimalbestfit_index_SPINEL,:,:};
            
            %saves primitve melt compositions
            BestFitPrimary_n(numdata,:) = CorrectedLiquids_zz(optimalbestfit_index,:);
            BestFitPrimary_Trace_n(numdata,:) = CorrectedLiquids_Trace_zz(optimalbestfit_index,:);
            
            correctedLiquids_optimal_SPINEL(numdata,:) = CorrectedLiquids_zz(optimalbestfit_index_SPINEL,:);
            
            
            %calculate phase transitions:
            Garnet2SpinelTransition(numdata,:) = (1.2-0.52.*BestFitPrimary_n(numdata, NaKNumCol)+2.11.*BestFitPrimary_n(numdata,MgNumCol))*10;
            Spinel2PlagTransition(numdata,:) = (0.823+1.139.*BestFitPrimary_n(numdata, NaKNumCol)+0.165.*BestFitPrimary_n(numdata,MgNumCol))*10;
            Garnet2SpinelTransition_SPIN(numdata,:) = (1.2-0.52.*correctedLiquids_optimal_SPINEL(numdata, NaKNumCol)+2.11.*correctedLiquids_optimal_SPINEL(numdata,MgNumCol))*10;
            Spinel2PlagTransition_SPIN(numdata,:) = (0.823+1.139.*correctedLiquids_optimal_SPINEL(numdata, NaKNumCol)+0.165.*correctedLiquids_optimal_SPINEL(numdata,MgNumCol))*10;
            
            %saves best fitting Mg#s at phase boundaries:
            %[fitness MgNumOPA MgNumOP];
            % [NRMSD	MG_CPX	MG_PLAG	FC extent	Mg# PM	Mg#OL NaK# PM	GarSpin	SpinPlag]
            SavedResults(numdata,:) = [optimalbestfit FCparameters_zz(optimalbestfit_index,:) CorrectedLiquids_zz(optimalbestfit_index,16) CorrectedLiquids_zz(optimalbestfit_index,[MgNumCol MgNumCol NaKNumCol]) Garnet2SpinelTransition(numdata,:) Spinel2PlagTransition(numdata,:)];
            SavedResults_SPINEL(numdata,:) = [optimalbestfit_SPINEL FCparameters_zz(optimalbestfit_index_SPINEL,:) CorrectedLiquids_zz(optimalbestfit_index_SPINEL,16) CorrectedLiquids_zz(optimalbestfit_index_SPINEL,[MgNumCol MgNumCol NaKNumCol]) Garnet2SpinelTransition_SPIN(numdata,:) Spinel2PlagTransition_SPIN(numdata,:)];
            
            %NRMSD-PM	MG_CPX	MG_PLAG	%OPC %OP %O	FC extent
            %Mg# PrimMelt	Mg# OLIV	NaK# PrimMelt
            SavedResults_v2(numdata,:)=[optimalbestfit FCparameters_zz(optimalbestfit_index,:) phaseChangeInfo_n(numdata,:)...
                CorrectedLiquids_zz(optimalbestfit_index,[MgNumCol MgNumCol NaKNumCol])];
            
            
            FCparameters_zz_NaN = FCparameters_zz(optimalbestfit_index,:);
            
            FCparameters_zz_NaN(FCparameters_zz_NaN==MinMgNum_Data) = NaN; % means the phase did not fractionate
            
            
            
            %
            %             SavedResults_v2_NaN(numdata,:)=[optimalbestfit FCparameters_zz_NaN phaseChangeInfo_n(numdata,:)...
            %                 CorrectedLiquids_zz(optimalbestfit_index,[MgNumCol MgNumCol NaKNumCol])];
            %
            
            
            SavedResults_v2_NaN(numdata,:)=[optimalbestfit FCparameters_zz_NaN phaseChangeInfo_n(numdata,:)...
                CorrectedLiquids_zz(optimalbestfit_index,[MgNumCol MgNumCol NaKNumCol])];
            
            
            Transitions_v2(numdata,:)=[Garnet2SpinelTransition(numdata,:) Spinel2PlagTransition(numdata,:)];
            
            %NRMSD-PM	MG_CPX	MG_PLAG	%OPC %OP %O	FC extent
            %Mg# PrimMelt	Mg# OLIV	NaK# PrimMelt
            SavedResults_v2_SPINEL(numdata,:)=[optimalbestfit_SPINEL FCparameters_zz(optimalbestfit_index_SPINEL,:) phaseChangeInfo_n(numdata,:) ...
                CorrectedLiquids_zz(optimalbestfit_index,[MgNumCol MgNumCol NaKNumCol]) Garnet2SpinelTransition(numdata,:) Spinel2PlagTransition(numdata,:)];
            
            
            
            %saves fit conditions [ NRMSD_GAR NRMSD_SPIN NRMSD_PLAG NRMSD_BEST BestFitFIELD]
            fits_optimal(numdata,:) = PrimitveMeltsFits_zz(optimalbestfit_index,:);
            fits_optimal_SPINEL(numdata,:) = PrimitveMeltsFits_zz(optimalbestfit_index_SPINEL,:);
            %saves fit conditions [ P_GAR T_GAR TMP_GAR  P_SPIN T_SPIN TMP_SPIN P_PLAG T_PLAG TMP_PLAG]
            fitsPTTmp_optimal(numdata,:) = PrimitveMelts_fitsPTTmp_zz(optimalbestfit_index,:);
            fitsPTTmp_optimal_SPINEL(numdata,:) = PrimitveMelts_fitsPTTmp_zz(optimalbestfit_index,4:6);
            
            % CorrectedLiquids_ff=
            % ['1-SiO2'	'2-TiO2' '3-Al2O3'	'4-Cr2O3'  '5-FeO'	'6-MnO'	'7-MgO'	'8-CaO'  '9-Na2O'	'10-K2O' '11-P2O5' '12-SUM 13-MG#MELT 14-MG#OLIV 15-NAK#MELT 16-FC_EXTENT]
            
            if fits_optimal(numdata,5) == 1
                AllFieldBestFit(numdata,:)=fitsPTTmp_optimal(numdata,1:3);
            else if fits_optimal(numdata,5)== 2
                    AllFieldBestFit(numdata,:)= fitsPTTmp_optimal(numdata,4:6);
                else if fits_optimal(numdata,5) == 3
                        AllFieldBestFit(numdata,:)= fitsPTTmp_optimal(numdata,7:9);
                    end
                end
            end
            
            
            OptimalPT(numdata,:) =[...
                BestFitP(numdata)...
                AllFieldBestFit(numdata,1:2)...
                ];
            
            
            %%%%% Saves information on all RFPs that come within error of
            %%%%% an MSP: (i.e., this is "Box g in Figure 4")
            
            GoodEggs_ReallyGood_ff_indicies = find(GoodEggs(:,5)<ErPrimary);
            GoodEggs_ReallyGood_All{numdata,:,:}=GoodEggs(GoodEggs_ReallyGood_ff_indicies,:);
            
            GoodEggs_ReallyGood_All_Average(numdata,:)=[mean(GoodEggs(GoodEggs_ReallyGood_ff_indicies,:),1) size(GoodEggs_ReallyGood_ff_indicies,1)];
            GoodEggs_ReallyGood_All_Stdev(numdata,:)=[std(GoodEggs(GoodEggs_ReallyGood_ff_indicies,:),0,1)];
            
            GoodEggs_All{numdata,:,:}=GoodEggs ;
            BagEggs_All{numdata,:,:}=BadEggs ;
            
            
            
            % Record how well each stability field fits the data:
            
            %[EachFieldMin_OPALM,EachFieldMin_OPALM_Index]=min(DistanceMatrix_NEW_ff_2,[],1);
            % https://www.mathworks.com/matlabcentral/answers/69517-find-first-value-below-a-minimum-in-a-vectorized-way
            A = EachFieldBestFit_zz_threshold(:,4:6)';
            thresh = ErPrimary;
            temp_mat = repmat(1:size(A,2),size(A,1),1);
            temp_mat = temp_mat .* (A<=thresh);
            temp_mat(temp_mat == 0) = NaN;
            [your_idx,your_idx2] = min(temp_mat,[],2);
            bestField_ff = your_idx';
            bestField_ff_index = your_idx2';
            
            
            [bestField_ff_ALL,bestField_ff_index_ALL]=min(EachFieldBestFit_zz_All(:,4:6),[],1);
            
            bestField_numdata_threshold{numdata,:,:}= [EachFieldBestFit_zz_threshold(bestField_ff_index(1),[2 3 4 7 10 13 ])
                EachFieldBestFit_zz_threshold(bestField_ff_index(2),[2 3 5 8 11 14])
                EachFieldBestFit_zz_threshold(bestField_ff_index(3),[2 3 6 9 12 15])];
            
            
            bestField_numdata_ALL{numdata,:,:}= [EachFieldBestFit_zz_All(bestField_ff_index_ALL(1),[2 3 4 7 10 13])
                EachFieldBestFit_zz_All(bestField_ff_index_ALL(2),[2 3 5 8 11 14])
                EachFieldBestFit_zz_All(bestField_ff_index_ALL(3),[2 3 6 9 12 15])];
            
            %                 AllCorrectedLiquids_temp_imported = XPTmeter_MASTER(BestFit_FCPath_n{numdata,:,:},TemperatureIndicies, CompIndicies, PressureIndicies, FittedVariablesIndicies, ModelVariablesIndicies,...
            %                     eval(sprintf('garnet_coefficients_%s',whichRegression)), eval(sprintf('spinel_coefficients_%s',whichRegression)), eval(sprintf('plagioclase_coefficients_%s',whichRegression)));
            
            
            
            % Variables with the entire best fit path for each basalt
            BestFit_FCPath_n;
            BestFit_FCPath_Tormey_n{numdata,:,:} = TormeyProjection(BestFit_FCPath_n{numdata,:,:},Elements,ElementIndicies4Target_Tormey);
            ReverseFCPath_Trace_n;
            
            clear DataRatios TraceElements
            TraceElements = ReverseFCPath_Trace_n{numdata,:,:};
            for gg = 1:size(ElementRatios,1)
                DataRatios(:,gg) = TraceElements(:,ElementRatiosIndex(gg,1)) ./TraceElements(:,ElementRatiosIndex(gg,2));
                DataRatios(:,21) = 2.*(TraceElements(:,EuCol)./TraceElements4Normalization_PUM(EuCol))./((TraceElements(:,SmCol)./TraceElements4Normalization_PUM(SmCol))+(TraceElements(:,GdCol)./TraceElements4Normalization_PUM(GdCol)));
                
            end
            ReverseFCPath_TraceRatios_n{numdata,:,:}=DataRatios;
            
            
            % Saves the best fit primary basalt for each erupted basalt
            BestFitPrimary_n;
            BestFitPrimaryTormey_n(numdata,:) = TormeyProjection(BestFitPrimary_n(numdata,:),Elements,ElementIndicies4Target_Tormey);
            BestFitPrimary_Trace_n;
            
            clear DataRatios TraceElements
            TraceElements = BestFitPrimary_Trace_n(numdata,:);
            for gg = 1:size(ElementRatios,1)
                DataRatios(:,gg) = TraceElements(:,ElementRatiosIndex(gg,1)) ./TraceElements(:,ElementRatiosIndex(gg,2));
                DataRatios(:,21) = 2.*(TraceElements(:,EuCol)./TraceElements4Normalization_PUM(EuCol))./((TraceElements(:,SmCol)./TraceElements4Normalization_PUM(SmCol))+(TraceElements(:,GdCol)./TraceElements4Normalization_PUM(GdCol)));
                
            end
            BestFitPrimary_TraceRatios_n(numdata,:)  = DataRatios;
            
            
            
            if optimalbestfit > ErPrimary && round(CorrectedLiquids_zz(optimalbestfit_index,MgNumCol),3)==MaxMgNumPrimary
                % Enter this loop if the best path has the same Mg# as the
                % max primary melt Mg# number but it farther from a MSP than ErPrimary, which means that the model
                % did not actually find a good fit
                GoodEggs_All{numdata,:,:}=GoodEggs ;
                BagEggs_All{numdata,:,:}=BadEggs ;
                
                SavedResults_v2(numdata,:)=NaN.*zeros(1,10);
                SavedResults_v2_NaN(numdata,:)=NaN.*zeros(1,10);
                fits_optimal(numdata,:)=NaN.*zeros(1,5);
                AllFieldBestFit(numdata,:)=NaN.*zeros(1,3);
                Transitions_v2(numdata,:)=NaN.*zeros(1,2);
                withinPhaseT(numdata,:)=NaN.*zeros(1,3);
                
                BestFitPrimary_n(numdata,:)=NaN.*zeros(1,30);
                BestFitPrimaryTormey_n(numdata,:)=NaN.*zeros(1,7);
                
                BestFitPrimary_Trace_n(numdata,:)=NaN.*zeros(1,45);
                BestFitPrimary_TraceRatios_n(numdata,:)=NaN.*zeros(1,22);
                BestFit_FCPath_Tormey_n(numdata,:)={NaN.*zeros(1,7)};
                BestFit_FCPath_n(numdata,:)={NaN.*zeros(1,30)};
                ReverseFCPath_Trace_n(numdata,:)={NaN.*zeros(1,45)};
                
                Garnet2SpinelTransition(numdata,:) = NaN; %(1.2-0.52.*BestFitPrimary_n(numdata, NaKNumCol)+2.11.*BestFitPrimary_n(numdata,MgNumCol))*10;
                Spinel2PlagTransition(numdata,:) = NaN; %(0.823+1.139.*BestFitPrimary_n(numdata, NaKNumCol)+0.165.*BestFitPrimary_n(numdata,MgNumCol))*10;
                
                FailMatrix{numdata}='Not Optimized';
            end
            
            
            if mod(numdata,50)==0
                eval(['save(''',[pwd '/' ImportPlotMasterFolder '/' subfolder '/' worksheetName2Save '.mat'],''')'])
            end
            
            
        end
        
    end
end
time=toc;

% This code just extracts more information for saving

%Check to make sure best-fit field is within pressure range expected for
% lherzolites

if exist('SavedResults')~=1
    disp('YIKES. All data failed. Try again...')
    SavedResults = NaN.*ones(numdata,5);
    FieldString = NaN.*ones(numdata,1);
    FieldPTString = NaN.*ones(numdata,1);
end

SaturationConditions=[SavedResults(:,2:3) SavedResults(:,5)] ;
fieldField = fits_optimal(:,5);
meltingP = AllFieldBestFit(:,1);

for d = 1:size(fieldField,1)
    if fieldField(d)==1
        if Transitions_v2(d,1)<meltingP(d)
            withinPhaseT(d,:) = [1 NaN meltingP(d)-Transitions_v2(d,1)];
            FieldPTString{d,1} = 'In Field';
            
        else
            withinPhaseT(d,:) = [0 meltingP(d)-Transitions_v2(d,1) NaN ];
            FieldPTString{d,1} = 'Out of Field';
        end
        FieldString{d,1} = 'Gar-lherz';
    end
    
    if fieldField(d)==2
        bothTransitions = abs([meltingP(d)-Transitions_v2(d,1) meltingP(d)-Transitions_v2(d,2)]);
        if meltingP(d)<Transitions_v2(d,1) && meltingP(d)>Transitions_v2(d,2)
            withinPhaseT(d,:) = [1 NaN min(bothTransitions)];
            FieldPTString{d,1} = 'In Field';
        else
            withinPhaseT(d,:) = [0 min(bothTransitions) NaN];
            FieldPTString{d,1} = 'Out of Field';
        end
        FieldString{d,1} = 'Spinel-lherz';
    end
    
    if fieldField(d)==3
        if meltingP(d)<Transitions_v2(d,2)
            withinPhaseT(d,:) = [1 NaN meltingP(d)-Transitions_v2(d,2)];
            FieldPTString{d,1} = 'In Field';
        else
            withinPhaseT(d,:) = [0 meltingP(d)-Transitions_v2(d,2) NaN];
            FieldPTString{d,1} = 'Out of Field';
        end
        FieldString{d,1} = 'Plag-lherz';
    end
    
    if isnan(fieldField(d))==1
        FieldPTString{d,1} = 'NaN';
        FieldString{d,1} = 'NaN';
    end
    
end

[ll]=find(withinPhaseT(:,1)==0 & abs(withinPhaseT(:,2))>2.5);
SavedResults_v2(ll,1)=0.1111;%SavedResults_v2(ll,1).*10;
SavedResults_v2_NaN(ll,1)=0.1111;%SavedResults_v2(ll,1).*10;

if any(ll)
    for ljmh = ll
        FailMatrix{ljmh}='!Outside Stability!';
    end
end


% Calculate equlibrum olivine composition at pressure of melting
Pb = AllFieldBestFit(:,1); % melting pressure in kbars
MgnumMelt = BestFitPrimary_n(:,MgNumCol);
KD_melting = 0.28 +.002.*Pb;%olivineKD_melting
MgnumOl_New = 1./(KD_melting.*(1./(MgnumMelt)-1)+1);

% Gather other useful parameters
MgErupt = data(:,MgNumCol);
FCP = BestFitP';
MgnumMelt = SavedResults_v2(:,8);
NaKMelt = SavedResults_v2(:,10);
PercentFC = SavedResults_v2(:,7);
MeltingTemp= AllFieldBestFit(:,2);
MantlePotentialTemp= AllFieldBestFit(:,3);
PM_error= SavedResults_v2(:,1);
KnownT =  runningMajors(:,1);
KnownP =  runningMajors(:,2);
KnownMgNum =  runningMajors(:,MgNumCol);

%Make OPALM_ERROR a string in percent
for z = 1:size(PM_error,1)
    PM_errorStr{z}=sprintf('%.1f%%',PM_error(z).*100);
end

%save nonunique info
izz=1;
for lhmg= 1:size(MantlePotentialTemp)
    GoodHere = GoodEggs_All{lhmg};
    ARedo_GoodEggs_All_indicies = find(GoodHere(:,5)<=0.015);
    GoodHere_reallyGood = GoodHere(ARedo_GoodEggs_All_indicies,:);
    GoodHere_reallyGood_Mean(izz,:) =[mean(GoodHere_reallyGood,1) size(ARedo_GoodEggs_All_indicies,1)];
    GoodHere_reallyGood_Stdev(izz,:) =[std(GoodHere_reallyGood,0,1)];
    izz=izz+1;
end
NonuniqueInfo = [GoodHere_reallyGood_Mean GoodHere_reallyGood_Stdev];

if strcmp(SaveFigs_OPALM_RMSD,'yes')
    savefigsPDF('',sprintf('Figs-OPALM-RMSD-%s',worksheetName2Save));
end


%% Make Matlab Tables!

OPALM_RMSD = SavedResults_v2_NaN(:,1);

TheBest = find(OPALM_RMSD<=0.015);

if size(subdata,2)==1
    RevPetResultsTable = table(round(subdata,0)',CombinedLabel,...
        AllFieldBestFit, {strcat(num2str(OPALM_RMSD*100,2),'%')'},FieldString,data(:,MgNumCol).*100,MgnumMelt.*100, MgnumOl_New.*100,...
        GoodHere_reallyGood_Mean(:,7), GoodHere_reallyGood_Stdev(:,7),...
        GoodHere_reallyGood_Mean(:,8), GoodHere_reallyGood_Stdev(:,8),...
        GoodHere_reallyGood_Mean(:,9), GoodHere_reallyGood_Stdev(:,9),...
        GoodHere_reallyGood_Mean(:,4).*100, GoodHere_reallyGood_Stdev(:,4).*100,...
        GoodHere_reallyGood_Mean(:,10), GoodHere_reallyGood_Stdev(:,10),...
        FieldPTString, withinPhaseT(:,2),FailMatrix',...
        Transitions_v2, withinPhaseT(:,3),...
        BestFitP',{strcat(num2str(OPAM_P_fitness_RMSD'*100,2),'%')},{strcat(num2str(PO_OPAM_fitness_RMSD'*100,2),'%')},...
        SavedResults_v2_NaN(:,2).*100, SavedResults_v2_NaN(:,3).*100,...
        {strcat(num2str(SavedResults_v2_NaN(:,4),3),'%')},{strcat(num2str(SavedResults_v2_NaN(:,5),3),'%')},{strcat(num2str(SavedResults_v2_NaN(:,6),3),'%')},...
        {strcat(num2str(SavedResults_v2_NaN(:,7),3),'%')},...
        BestFitPrimary_n(:,3:end), BestFitPrimary_Trace_n,BestFitPrimary_TraceRatios_n);
else
    RevPetResultsTable = table(round(subdata,0)',CombinedLabel,...
        AllFieldBestFit, strcat(num2str(OPALM_RMSD*100,2),'%'),FieldString,data(:,MgNumCol).*100,MgnumMelt.*100, MgnumOl_New.*100,...
        GoodHere_reallyGood_Mean(:,7), GoodHere_reallyGood_Stdev(:,7),...
        GoodHere_reallyGood_Mean(:,8), GoodHere_reallyGood_Stdev(:,8),...
        GoodHere_reallyGood_Mean(:,9), GoodHere_reallyGood_Stdev(:,9),...
        GoodHere_reallyGood_Mean(:,4).*100, GoodHere_reallyGood_Stdev(:,4).*100,...
        GoodHere_reallyGood_Mean(:,10), GoodHere_reallyGood_Stdev(:,10),...
        FieldPTString, withinPhaseT(:,2),FailMatrix',...
        Transitions_v2, withinPhaseT(:,3),...
        BestFitP',strcat(num2str(OPAM_P_fitness_RMSD'*100,2),'%'),strcat(num2str(PO_OPAM_fitness_RMSD'*100,2),'%'),...
        SavedResults_v2_NaN(:,2).*100, SavedResults_v2_NaN(:,3).*100,...
        strcat(num2str(SavedResults_v2_NaN(:,4),3),'%'),strcat(num2str(SavedResults_v2_NaN(:,5),3),'%'),strcat(num2str(SavedResults_v2_NaN(:,6),3),'%'),...
        strcat(num2str(SavedResults_v2_NaN(:,7),3),'%'),...
        BestFitPrimary_n(:,3:end), BestFitPrimary_Trace_n,BestFitPrimary_TraceRatios_n);
end


% Save significant digits
numericVars = find(varfun(@isnumeric,RevPetResultsTable,'output','uniform'));
RevPetResultsTable(:,numericVars) = varfun(@(var) round(var, 4, 'significant'), RevPetResultsTable(:,numericVars));


%TabName','DataName','Region_Label'

RevPetResultsTable.Properties.VariableNames ={'PointNum','Region_Label','AllFieldBestFit','OPALM_RMSD','StabilityField','MgNum_Erupted', 'MgNum_Primary','MgNum_Olivine',...
    'Avg_P_kbar','Stdev_P_kbar','Average_T_C','Stdev_T_C','Average_TP_C','Stdev_TP_C',...
    'Avg_MgNum','Stdev_MgNum','Average_percentFC','Stdev_percentFC','FieldPTString','DistanceFromField','FailureReason','GarSpPlagTransitions','DistanceToNearestBoundary'...
    'PressureFC','OPAM_RMSD_Erupted','OPM_RMSD_Erupted','MG_CPX','MG_PLAG','OPA_FC','OP_FC','O_FC','totalFC','PrimaryMajor','PrimaryTrace','PrimaryRatios'};

NewElements = {'SiO2_Primary'	'TiO2_Primary'	  'Al2O3_Primary' 'Cr2O3_Primary'	  'FeO_Primary'	  'MnO_Primary'	  'MgO_Primary'	  'CaO_Primary' 'Na2O_Primary'   'K2O_Primary' 'P2O5_Primary' 'NiO_Primary' 'H2O_Prim' 'total_Prim'...
    'MgNum_Prim' 'NaKNum_Primary' 'CaOAl2O3wt_Primary' 'CaOAl2O3mol_Primary' 'Na2OFeO_Primary' 'K2OTiO2_Primary' 'OneMinusMgNum_Primary' 'Qtz_Primary' 'Plag_Primary' 'Oliv_Primary' 'Cpx_Primary' 'Ox_Primary' 'Or_Primary' 'Ap_Primary'};



RevPetResultsTable = splitvars(RevPetResultsTable,'AllFieldBestFit','NewVariableNames',{'Melting_P_kbar' 'Melting_T_C','apparentT_P'});
RevPetResultsTable = splitvars(RevPetResultsTable,'PrimaryMajor','NewVariableNames',NewElements);
RevPetResultsTable = splitvars(RevPetResultsTable,'PrimaryTrace','NewVariableNames', strcat(targetStrings_Trace,'_Primary'));
RevPetResultsTable = splitvars(RevPetResultsTable,'PrimaryRatios','NewVariableNames', strcat(NewRatioLabels,'_Primary'));

%%%%%%%%%%%%%%%%%%%%

% this is just here to deal with 1 versus 1+ data points in the table,
% which have to be formatted differently
RevPetResultsTable_ResultsOnly = RevPetResultsTable;

if size(subdata,2)==1
    RevPetResultsTable_ResultsOnly.OPALM_RMSD =  str2double(regexprep(RevPetResultsTable_ResultsOnly.OPALM_RMSD{:}','%',''));
    RevPetResultsTable_ResultsOnly.totalFC =  str2double(regexprep(RevPetResultsTable_ResultsOnly.totalFC{:},'%',''));
    RevPetResultsTable_ResultsOnly.OPA_FC =  str2double(regexprep(RevPetResultsTable_ResultsOnly.OPA_FC{:},'%',''));
    RevPetResultsTable_ResultsOnly.OP_FC =  str2double(regexprep(RevPetResultsTable_ResultsOnly.OP_FC{:},'%',''));
    RevPetResultsTable_ResultsOnly.O_FC =  str2double(regexprep(RevPetResultsTable_ResultsOnly.O_FC{:},'%',''));
    RevPetResultsTable_ResultsOnly.OPAM_RMSD_Erupted =  str2double(regexprep(RevPetResultsTable_ResultsOnly.OPAM_RMSD_Erupted{:},'%',''));
    RevPetResultsTable_ResultsOnly.OPM_RMSD_Erupted =  str2double(regexprep(RevPetResultsTable_ResultsOnly.OPM_RMSD_Erupted{:},'%',''));
else
    RevPetResultsTable_ResultsOnly.OPALM_RMSD = str2double(cellfun(@(x) regexprep(x,'%',''),cellstr(RevPetResultsTable_ResultsOnly.OPALM_RMSD),'uni',0));
    RevPetResultsTable_ResultsOnly.totalFC = str2double(cellfun(@(x) regexprep(x,'%',''),cellstr(RevPetResultsTable_ResultsOnly.totalFC),'uni',0));
    RevPetResultsTable_ResultsOnly.OPA_FC = str2double(cellfun(@(x) regexprep(x,'%',''),cellstr(RevPetResultsTable_ResultsOnly.OPA_FC),'uni',0));
    RevPetResultsTable_ResultsOnly.OP_FC = str2double(cellfun(@(x) regexprep(x,'%',''),cellstr(RevPetResultsTable_ResultsOnly.OP_FC),'uni',0));
    RevPetResultsTable_ResultsOnly.O_FC = str2double(cellfun(@(x) regexprep(x,'%',''),cellstr(RevPetResultsTable_ResultsOnly.O_FC),'uni',0));
    RevPetResultsTable_ResultsOnly.OPAM_RMSD_Erupted = str2double(cellfun(@(x) regexprep(x,'%',''),cellstr(RevPetResultsTable_ResultsOnly.OPAM_RMSD_Erupted),'uni',0));
    RevPetResultsTable_ResultsOnly.OPM_RMSD_Erupted = str2double(cellfun(@(x) regexprep(x,'%',''),cellstr(RevPetResultsTable_ResultsOnly.OPM_RMSD_Erupted),'uni',0));
end
RevPetResultsTable=[NewTable(:,6:9) RevPetResultsTable(:,1:34),NewTable(:,[1:4,10:20]),RevPetResultsTable(:,35:end), NewTable(:,21:end)];
%% write the table to a spreadsheet

filename = sprintf('%s.xlsx',worksheetName2Save);
writetable(RevPetResultsTable,[pwd '/' ImportPlotMasterFolder '/' subfolder '/' filename],'Sheet','DetailedResults','UseExcel',1)

NumberRun = size(runningMajors,1);
NumberWellFit= size(TheBest,1);

AvgStdMinMaxResults=[mean(RevPetResultsTable{:,7}) std(RevPetResultsTable{:,7}) min(RevPetResultsTable{:,7}) max(RevPetResultsTable{:,7});...
    mean(RevPetResultsTable{:,8}) std(RevPetResultsTable{:,8}) min(RevPetResultsTable{:,8}) max(RevPetResultsTable{:,8});...
    mean(RevPetResultsTable{:,9}) std(RevPetResultsTable{:,9}) min(RevPetResultsTable{:,9}) max(RevPetResultsTable{:,9});...
    mean(OPALM_RMSD) std(OPALM_RMSD) min(OPALM_RMSD) max(OPALM_RMSD)]';
AvgStdMinMaxResults=[AvgStdMinMaxResults; NumberRun.*ones(1,4)];



AvgStdMinMaxBestResults=[mean(RevPetResultsTable{TheBest,7}) std(RevPetResultsTable{TheBest,7}) min(RevPetResultsTable{TheBest,7}) max(RevPetResultsTable{TheBest,7});...
    mean(RevPetResultsTable{TheBest,8}) std(RevPetResultsTable{TheBest,8}) min(RevPetResultsTable{TheBest,8}) max(RevPetResultsTable{TheBest,8});...
    mean(RevPetResultsTable{TheBest,9}) std(RevPetResultsTable{TheBest,9}) min(RevPetResultsTable{TheBest,9}) max(RevPetResultsTable{TheBest,9});...
    mean(OPALM_RMSD(TheBest)) std(OPALM_RMSD(TheBest)) min(OPALM_RMSD(TheBest)) max(OPALM_RMSD(TheBest))];

if size(AvgStdMinMaxBestResults,2)~=4
    AvgStdMinMaxBestResults=NaN.*ones(4,4);
end

AvgStdMinMaxBestResults=[AvgStdMinMaxBestResults; NumberWellFit.*ones(1,4)];


ROWS = {'Mean','Standard Deviation','Min','Max','Number'}';
AveragedResults = table(ROWS,AvgStdMinMaxResults,AvgStdMinMaxBestResults);


AveragedResults = splitvars(AveragedResults,'AvgStdMinMaxResults','NewVariableNames', {'All_MeltingP_kbar','All_MeltingT_C','All_MeltingApparentTP_C','All_OPALM_RMSD'});
AveragedResults = splitvars(AveragedResults,'AvgStdMinMaxBestResults','NewVariableNames', {'Best_MeltingP_kbar','Best_MeltingT_C','Best_MeltingApparentTP_C','Best_OPALM_RMSD'});
writetable(AveragedResults,[pwd '/' ImportPlotMasterFolder '/' subfolder '/' filename],'Sheet','AverageResults','UseExcel',1)

writetable(NewTable,[pwd '/' ImportPlotMasterFolder '/' subfolder '/' filename],'Sheet','ReformattedData','UseExcel',1)



inputvalues = table({whichRegression},MinMgNumPrimary,MaxMgNumPrimary, MinMgNumOP ,MaxMgNumOP, MinMgNumOPA, MaxMgNumOPA, ErOPC, ErOP, ErPrimary,time,NumberRun,NumberWellFit);
inputvalues.Properties.VariableNames{'time'} = 'TimeinSec';
writetable(inputvalues,[pwd '/' ImportPlotMasterFolder '/' subfolder '/' filename],'Sheet','Settings','UseExcel',1)


% %MAYBE
% 'Gar fit','Sp fit','Plag fit',
% fits_optimal(:,1:4)...
% 'Gar fit','Sp fit','Plag fit','Best fit',...


eval(['save(''',[pwd '/' ImportPlotMasterFolder '/' subfolder '/' worksheetName2Save '.mat'],''')'])

'Name of Folder Files are Saved in:'
subfolder
'Name of Saved Matlab Workspace and EXCEL spreadsheet:'
worksheetName2Save


toc
load train;
sound(y,1/2*Fs)
